## Script simulation PA note
## Last modified: 2014-01-06

## Start
rm(list=ls())
library(foreach); library(doSNOW); require(foreign)
registerDoSNOW(makeCluster(3, type = "SOCK"))

## Set seed
set.seed(77137248) # a dollar bill in a friend's pocket...

## Sims fcn in parallel...
# nsim =  # of simulations; Ns = # of obs; Ts = Time periods; 
# rhos = vector of rhos and phis
# burn =  data to burn; gamma = parameter in the dgp for 
# Y_t = \phi Y_t-1 + gamma*X_t-1 + epsilon_t
# dir = directory to save stata files
fsim <- function(nsim, Ns, Ts, rhos, iFEx, iFEy, burn=100, gamma=0, dir) {
  # write result in a stata file
  write_stata <- function (res, dir) {
    withFEx <- function(x) ifelse(x, 'withFEx', 'woFEx')
    for (j in 1:length(res)) {
      filename = paste(withFEx(res[[j]][1,'FEx']),
                       'rho', res[[j]][1,'rho'],
                       'phi', res[[j]][1,'phi'],
                       'gamma', res[[j]][1,'gamma'],
                       'N', res[[j]][1,'n'],
                       'T', paste(res[[j]][1,'t'],
                                 '.dta', sep=''),
                       sep='_')
      write.dta(res[[j]][,c("rho_hat","b_hat", "Rej_b_diff0",
                            "se_rho_hat", "se_b_hat")], 
                paste(dir, filename, sep='/'))
    }
  }
  # DGP
  build_tscs_ar1 <- function (rho, Nind, Ti, FEy, FEx, phi, burn, gamma) {
    res <- foreach (i = 1:Nind, .combine=rbind) %dopar% {
      # DGP for Y and X
      nu = rnorm(Ti+burn, mean=0, sd=1)
      epsilon = rnorm(Ti+burn, mean=0, sd=1)
      Y = rep(0, Ti+burn)    
      X = rep(0, Ti+burn)
      X[1] = nu[1]
      Y[1] = epsilon[1]
      for (j in 2:(Ti+burn)) {
        Y[j] = FEy[i] + rho*Y[j-1] + gamma*X[j-1] + epsilon[j]
        X[j] = FEx[i] + rho*X[j-1] + phi*Y[j-1] + nu[j]
      }
      # Burn the initialization
      X = X[-c(1:burn-1)]
      Y = Y[-c(1:burn-1)]
      # Compute lags
      X_lag1 = c(NA,X[-length(X)])
      Y_lag1 = c(NA,Y[-length(Y)])
      aux=data.frame(Ind=rep(i,Ti+1), Time=0:Ti, Y=Y, 
                     Y_lag1 = Y_lag1, X=X, X_lag1=X_lag1)
      aux
    }
    res <- data.frame(res)
    return(res)
  }
  # Beginning
  regs <- list()
  for (i in rhos) {
    cat('\n', 'T = ', Ts)
    cat('\n', 'rho = ', i[1], ' phi = ', i[2], '\n', 'N = ')
    for (k in Ns) {
      cat(k, '\t')
      FEx = rep(0,k)
      FEy = rep(0,k)
      r <- foreach (l = 1:nsim, .combine=rbind, .packages='foreach') %dopar% {
        if (iFEx) FEx = rnorm(k, mean=0, sd=1)
        if (iFEy) FEy = rnorm(k, mean=0, sd=1)
        simdata <- build_tscs_ar1(i[1], k, Ts, FEy, FEx, phi=i[2], burn, gamma)
        if (iFEx) model <- lm(Y~Y_lag1+X_lag1+as.factor(Ind), data=simdata)
        else model <- lm(Y~Y_lag1+X_lag1, data=simdata)
        rm(simdata)
        c(i[1], i[2], k, Ts, model$coef[2:3],
          as.numeric(abs(summary(model)$coef[3,3])>1.96),
          summary(model)$coef[2,2], summary(model)$coef[3,2], 
          iFEx, iFEy, gamma)
      }
      r <- data.frame(r)
      names(r) <- c('rho', 'phi', 'n', 't', 'rho_hat', 'b_hat', 'Rej_b_diff0', 
                    'se_rho_hat', 'se_b_hat', 'FEx', 'FEy', 'gamma')
      row.names(r) <- NULL
      r <- list(r)
      regs <- append(regs, r)
    }
  }
  cat('\n','Finish... \n')
  write_stata(regs, dir)
}

# Initialization
nsim = 1000 # Number of sims
Ns = c(20, 40, 80, 200) # N

cat('\n', 'Starting simulations...')

## First sim: like todd...
rhos = list(c(.3,.3))
system.time(fsim(nsim, Ns, Ts=40, rhos, iFEx = T, iFEy = T, burn=100, gamma=.3, dir='~/Desktop'))[3]
system.time(fsim(nsim, Ns, Ts=400, rhos, iFEx = T, iFEy = T, burn=100, gamma=.3, dir='~/Desktop'))[3]


cat('\n', 'End of script...')
## End of script
